home *** CD-ROM | disk | FTP | other *** search
/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc A) / Acorn User China CD-ROM (UK) (Disc A).bin / HENSA / MATHS / JACAL.ARC / jacal_vect < prev    next >
Encoding:
Text File  |  1993-08-27  |  6.7 KB  |  266 lines
  1. ;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
  2. ;;; Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
  3. ;;; See the file "COPYING" for terms applying to this program.
  4.  
  5. (define (bunch:norm x)
  6.   (if (null? (cdr x))
  7.       (car x)
  8.     x))
  9.  
  10. (define (copymatrix m)
  11.   (if (matrix? m)
  12.       (map copy-list m)
  13.       (math:error 'Not-a-matrix:- m)))
  14.  
  15. (define (row? a)
  16.   (and (bunch? a)
  17.        (notevery bunch? a)))
  18.  
  19. (define (matrix? a)
  20.   (and (bunch? a)
  21.        (not (null? a))
  22.        (bunch? (car a))
  23.        (let ((len (length (car a))))
  24.      (and (every (lambda (r) (and (bunch? r) (= (length r) len)))
  25.              (cdr a))
  26.           len))))
  27.  
  28. (define (matrix . lst)
  29.   (if (matrix? lst) lst (matrix lst)))
  30.  
  31. (define (butnth n lst)
  32.   (cond ((null? lst) (math:error 'Coordinate-out-of-range:- n " " lst))
  33.     ((zero? n) (cdr lst))
  34.     (else (cons (car lst)
  35.             (butnth (+ -1 n) (cdr lst))))))
  36.  
  37. (define (mtrx:minor m i j)
  38.   (butnth (+ -1  i)
  39.       (map (lambda (x) (butnth (+ -1 j) x))
  40.            m)))
  41.  
  42. (define (transpose m)
  43.   (cond ((matrix? m)
  44.      (bunch:norm (apply map list m)))
  45.     ((bunch? m) (map list m))
  46.     (else m)))
  47.  
  48. (define (cart-prod choices)
  49.   (if (null? choices) '(())
  50.       (let* ((choice (car choices)))
  51.     (apply append
  52.            (map (lambda (tuple)
  53.               (map (lambda (elt)
  54.                  (cons elt tuple))
  55.                choice))
  56.             (cart-prod (cdr choices)))))))
  57.  
  58. (define (mtrx:square? m)
  59.   (if (not (matrix? m))
  60.       #f
  61.     (let ((rows (length m)))
  62.       (every (lambda (r) (= (length r) rows)) m))))
  63.  
  64. (define (mtrx:genmatrix fun i2 j2 i1 j1)
  65.   (do ((i i2 (+ -1 i))
  66.        (mat '()
  67.         (cons
  68.          (do ((j j2 (+ -1 j))
  69.           (row '() (cons (app* fun i j) row)))
  70.          ((< j j1) row))
  71.          mat)))
  72.       ((< i i1)
  73.        (bunch:norm mat))))
  74.  
  75. ;;;Create a matrix of indeterminates
  76.  
  77. (define (nmatrix var . rst)
  78.   (apply genmatrix
  79.      (lambda subs (apply rapply var subs))
  80.      rst))
  81.  
  82. (define (dotproduct a b)
  83.   (let ((cols (matrix? a))
  84.     (colbs (matrix? b)))
  85.     (cond ((or (and cols (bunch? b)) (and colbs (bunch? a)))
  86.        (if (not (and (eqv? cols colbs) (eqv? (length a) (length b))))
  87.            (math:error 'dotproduct-of-unequal-size-matrices:- a
  88.                " " b))
  89.        (reduce (lambda (x y) (app* $1+$2 x y))
  90.            (reduce (lambda (x y) (app* $1+$2 x y))
  91.                (map (lambda (x y) (app* $1*$2 x y))
  92.                 a b))))
  93.       ((and (row? a) (row? b))
  94.        (reduce (lambda (x y) (app* $1+$2 x y))
  95.            (map (lambda (x y) (app* $1*$2 x y))
  96.             a b)))
  97.       ((bunch? a) (map (lambda (x) (app* $1*$2 x b))
  98.                a))
  99.       ((bunch? b) (map (lambda (x) (app* $1*$2 a x))
  100.                b))
  101.       (else (app* $1*$2 a b)))))
  102.  
  103. (define (ncmult a b)
  104.   (let ((cols (matrix? a))
  105.     (colbs (matrix? b)))
  106.     (cond ((or (and cols (bunch? b)) (and colbs (bunch? a)))
  107.        (if (not colbs) (set! b (list b)))
  108.        (if (not (= (or cols (length a)) (length b)))
  109.            (math:error 'matrix-product-of-unequal-size-matrices:- a
  110.                " " b))
  111.        (or cols (set! a (matrix a)))
  112.        (bunch:norm
  113.         (map (lambda (arow)
  114.            (apply map
  115.               (lambda bcol
  116.                 (dotproduct bcol arow))
  117.               b))
  118.          a)))
  119.       (else (deferop _ncmult a b)))))
  120.  
  121. (define (mtrx:expt a pow)
  122.   (cond ((zero? pow) 1)
  123.     ((one? pow) a)
  124.     ((negative? pow) (mtrx:expt (mtrx:inverse a) (- pow)))
  125.            (else (ipow-by-squaring a (+ -1 pow) a ncmult))))
  126.  
  127. (define (crossproduct a b)
  128.   (if (and (row? a) (row? b))
  129.       (cond ((not (= (length a) (length b)))
  130.          (math:error 'Crossproduct 'unequal-length-vectors a b))
  131.         ((= 2 (length a))
  132.          (app* $1-$2
  133.            (app* $1*$2 (car a) (cadr b))
  134.            (app* $1*$2 (car b) (cadr a))))
  135.         ((= 3 (length a))
  136.          (list
  137.           (app* $1-$2
  138.             (app* $1*$2 (cadr a) (caddr b))
  139.             (app* $1*$2 (cadr b) (caddr a)))
  140.           (app* $1-$2
  141.             (app* $1*$2 (caddr a) (car b))
  142.             (app* $1*$2 (caddr b) (car a)))
  143.           (app* $1-$2
  144.             (app* $1*$2 (car a) (cadr b))
  145.             (app* $1*$2 (car b) (cadr a)))))
  146.         (else (math:error 'Crossproduct 'funny-length-vectors a b)))
  147.       (dotproduct a b)))
  148.  
  149. (define (mtrx:scalarmatrix n x)
  150.   (do ((i 1 (+ 1 i))
  151.        (rows (list (nconc (make-list (+ -1 n) 0) (list x)))
  152.          (cons (append (cdar rows) (list 0)) rows)))
  153.       ((>= i n) rows)))
  154.  
  155. (define (mtrx:diagmatrix diag)
  156.   (set! diag (reverse diag))
  157.   (do ((i (+ -1 (length diag)) (+ -1 i))
  158.        (j 0 (+ 1 j))
  159.        (diag diag (cdr diag))
  160.        (rows '()
  161.          (cons (nconc (make-list i 0) (list (car diag)) (make-list j 0))
  162.            rows)))
  163.       ((null? diag) rows)))
  164.  
  165. (define (determinant m)
  166.   (if (not (mtrx:square? m)) (math:error 'determinant-of-non-square-matrix m))
  167.   (letrec
  168.       ((x11 1)
  169.        (ds (lambda (m)
  170.          (cond ((null? m) 0)
  171.            ((expr:0? (caar m))
  172.             (set! x11 (- x11))
  173.             (cons (cdar m) (ds (cdr m))))
  174.            (else
  175.             (set! x11
  176.               (if (negative? x11) (app* _-$1 (caar m)) (caar m)))
  177.             (map (lambda (ro)
  178.                (if (expr:0? (car ro))
  179.                    (cdr ro)
  180.                    (app* $1-$2*$3
  181.                      (cdr ro)
  182.                      (cdar m)
  183.                      (app* $1/$2 (car ro) (caar m)))))
  184.              (cdr m)))))))
  185.     (cond ((null? (cdr m)) (caar m))
  186.       ((null? (cddr m)) (crossproduct (car m) (cadr m)))
  187.       (else (let ((subdet (determinant (ds m))))
  188.           (app* $1*$2 x11 subdet))))))
  189.  
  190. (define (charpoly m var)
  191.   (determinant (app* $1-$2
  192.              m
  193.              (mtrx:scalarmatrix (length m) var))))
  194.  
  195. (define (cofactor m i j)
  196.   (let ((det (determinant (mtrx:minor m i j))))
  197.     (if (odd? (+ i j))
  198.     (app* _-$1 det)
  199.     det)))
  200.  
  201. (define (coefmatrix eqns vars)
  202.   (map
  203.    (lambda (eq)
  204.      (map (lambda (var)
  205.         (let ((c (poly:coeff eq var 1)))
  206.           (set! eq (poly:coeff eq var 0))
  207.           c))
  208.       vars))
  209.    eqns))
  210.  
  211. (define (augcoefmatrix eqns vars)
  212.   (map
  213.    (lambda (eq)
  214.      (nconc
  215.       (map (lambda (var)
  216.          (let ((value (poly:coeff eq var 1)))
  217.            (set! eq (poly:coeff eq var 0))
  218.            value))
  219.        vars)
  220.       (list eq)))
  221.    eqns))
  222.  
  223. ;;; This algorithm taken from:
  224. ;;; Knuth, D. E.,
  225. ;;; The Art Of Computer Programming, Vol. 2: Seminumerical Algorithms,
  226. ;;; Addison Wesley, Reading, MA 1969.
  227. ;;; pp 425-426
  228. (define (rank m)
  229.   (let* ((cols (matrix? m))
  230.      (n (length m))
  231.      (c (make-list n -1))
  232.      (r 0))
  233.     (cond ((null? cols) (set! m (list m))
  234.             (set! n 1))
  235.       ((> cols n) (set! m (copymatrix m)))
  236.       (else (set! m (transpose m))
  237.         (set! n cols)))
  238.     (do ((k 0 (+ 1 k)))
  239.     ((>= k n) (- n r))
  240.       (let* ((l #f)
  241.          (j
  242.           (do ((j 0 (+ 1 j)))
  243.           ((or l (>= j n)) l)
  244.         (if (and (< (list-ref c j) 0)
  245.              (not (poly:0? (list-ref (list-ref m k) j))))
  246.             (set! l j))))
  247.          (rowj (and j (list-ref m j))))
  248.     (cond (j
  249.            (set! rowj (app* _-$1/$2
  250.                 rowj
  251.                 (list-ref rowj k)))
  252.            (set-car! (list-tail m j) rowj)
  253.            (do ((restrows m (cdr restrows)))
  254.            ((null? restrows))
  255.          (if (not (eq? (car restrows) rowj))
  256.              (set-car! restrows
  257.                    (app* $1*$2+$3
  258.                      rowj
  259.                      (list-ref (car restrows) k)
  260.                      (car restrows)))))
  261.            (set-car! (list-tail c j) k))
  262.           (else (set! r (+ 1 r))))))))
  263.  
  264. ;;;    Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
  265.  
  266.